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We discuss the structure and asymptotic long-time properties ol coupled equations for the mo- 
ments of a Brownian particle's momentum {p"'{t)) derived microscopically beyond the lowest ap- 
proximation in the weak coupling parameter A. Generalized fluctuation-dissipation relations are 
derived and shown to ensure convergence to thermal equilibrium to any order in A. 
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I. INTRODUCTION 



Many popular results of non-equilibrium statistical me- 
chanics, such as exponential decay of correlations, hold 
only on a time scale much longer than the characteristic 
relaxation time Tf, for a thermal bath and are thus ap- 
plicable only for sufficiently slow systems with relaxation 
time Ts ^ Tb- Such results may be obtained microscopi- 
cally in the lowest order of the perturbation theory with 
the ratio A ^ t^/ts as a small parameter and, addition- 
ally, with a coarse-grained time resolution much larger 
than Tfo. This approximation is referred to as the weak 
coupling limit and can be concisely formulated as a com- 
bination of three conditions A — ^ 0, t — > oo, with At 
finite. 

There are, of course, many situations of physical in- 
terest when the weak coupling approximation does not 
suffice [ll-Hl- Van Kampen developed a highly success- 
ful theory which allows us to take into account terms of 
higher orders in A, but still using as a prerequisite the 
time coarse graining in a form of the assumption that 
the system interacts with the bath via 'instantaneous' 
binary collisions The relaxing of this rather artifi- 
cial limit leads, in general, to non-Markovian master or 
Langevin equations which are more difficult to handle 
than their Markovian counterparts. Beyond the weak 
coupling limit, these equations have a rather compli- 
cated structure Q , and not much is known about their 
properties. In particular, while van Kampen's theory is 
consistent with equilibrium statistical mechanics [6], the 
relaxation to Maxwell-Boltzmann equilibrium within a 
more general approach, which takes into account multi- 
ple collisions and non-Markovian effects, is not entirel 
obvious and was questioned in several studies (see 0, 
and references therein). 

One purpose of this paper is to put the equations of 
stochastic dynamics into a form convenient for the eval- 
uation of stationary solutions to any order in the weak 
coupling parameter A. Most previous works in this di- 
rection seek to generalize the Fokker-Planck equation for 
the distribution function f{a,t) for a targeted dynami- 
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cal variable a. One difficulty with this approach is that 
beyond the lowest order in A generalized Fokker-Planck 
equations involve derivatives 9" //9a" of orders higher 
than two and do not guarantee positive definiteness of 
the solution. Also, and perhaps more importantly, within 
this approach it is not clear how to, in a systematic per- 
turbative way, take into account non-Markovian effects. 
The same perturbation technique which justifies results 
in the weak coupling limit may not work well in higher 
orders in A, leading to stationary distributions for the 
system with an equilibrium temperature different than 
that of the bath Q . Neither real @ nor numerical 
experiments suggest such a possibility. Some other trou- 
bles with non-Markovian Fokker-Planck equations were 
reported in [IT'] . In this paper we show that some of these 
difficulties can be avoided if one works with microscopi- 
cally derived Langevin equations for the powers a" of a 
targeted variable. These equations can be readily used to 
derive equations for the moments (a"(t)) which are lin- 
ear and not difficult to work with even in non-Markovian 
form, at least as far as stationary solutions are concerned. 

Although the discussion can be carried on a very gen- 
eral level, we choose to consider, for the sake of better 
visualization, the archetype example of a Brownian par- 
ticle of mass M immersed in a infinitely large bath at 
temperature T = 1/(3 and composed of molecules with 
mass m <C M. We will employ the method by Albers 
et al. (l2j to derive equations of motion for the moments 
(p"(t)) of the Brownian particle's (scaled) momentum. It 
is shown that to any order in the weak coupling parame- 
ter A = (to/M)^/^ the moments relax to the equilibrium 
values prescribed by equilibrium statistical mechanics. 
Convergence to thermal equilibrium is guided by gener- 
alized fluctuation-dissipation relations, which can hardly 
be derived by any other method but microscopically. On 
the other hand, thermahzation is found to be insensitive 
to particular relations between involved microscopic cor- 
relation functions. This leads to the optimistic conclusion 
that a consistent theory need not be totally microscopic. 



II. EXACT EQUATIONS 

Consider a system of N bath molecules of mass m in- 
teracting with each other and with a Brownian particle 
of mass M via short range potential U. The Hamiltonian 
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of the system is _ff = P^/2M + Ho{x), where P is the 
momentum of the particle, and Hq{x) is the Hamihonian 
of bath molecules in the field of the particle fixed at the 
position X. Introducing as usual the scaled momentum 
of the particle p = \P with A = (m/M)^/^, the Liouville 
operator of the systems can be written as 



L = Lo + XLi 



(1) 



where the operator Li acts on the particle's variables 
only 



_ p d d 
m ox op 



(2) 



while Lq governs the dynamics of the bath in the field of 
the fixed particle, 



N 



^ |_ m dxi Opi 



(3) 



and thus satisfies the relation LqHq ~ 0. In these 
equations Xi, Pi are coordinates and momenta of bath 
molecules, Fi = —dU/dxi and F = —dU/dx are the 
forces on the ith molecules and on the particle, respec- 
tively. 

We wish to decompose the exact equation 



dt 



(4) 



where p" — p"(0), into a form convenient to derive the 
Langevin equation for p" using an expansion in the small 
parameter A. Using the operator identity 

"'0 

with ^ = L, B = -VL, and Q = 1 - 75 one gets 



"'0 



(6) 



The operator V is convenient to chose to be a projector 
operator iV'^ = V) that averages over initial values of 
bath variables 



VA={A)= j p Adxi...dxNdpi...dpN 



with the canonical distribution 

1 



P = — e 
Z 



(7) 



(8) 



Such defined projection operator V and its complement 
2 = 1 — 7^ satisfy the relations 

VL^ = 0, QLo - Lo. (9) 

With ^ and the right-hand side of @ can be written 



Then (|4]) takes the desirable pre-Langevin form 

= \KAt) + \^ /*dre^(*-^)7'Lii^„(T), (10) 
at Jo 



where 



(11) 



plays the role of a "random" (rapidly fluctuating) force. 

Notice that the above expression for Kn{t) can be al- 
ternatively written with an additional factor 2 = 1 — 7': 

Knit) = A-ie2^*QLp". (12) 

This is because 

T'Lp" = XVLip"^ = Xnp'^-^VF = Xnp'^-^F) = 0. 



The form p2|) makes it obvious that the fluctuating term 
K{t) is zero centered 



{K,,it))^VKjt)^VQ^O. 



(13) 



Another useful identity involving Kn{t), which can be 
readily proved by integrating by parts, reads 



-g^Knit)) = 'P-Q^K^t) = -PiFKnit)}. (14) 



Equipped with these relations, one eventually write Eq. 
in the form 



dt 



p"(t) = XKnit) 



(15) 



A' / dT 



/ \dp m 

with the fluctuating force 

Kn{t)=ne^^^ Fp^-^. (16) 

The procedure outlined above is generic and can be 
easily generalized to derive exact Langevin-like equation 
of motion for an arbitrary dynamical variable or distri- 
bution function [12i] . 

III. LOWEST-ORDER PERTURBATION 

As it is, the exact equation (|15p is of little help because 
it contains the variables of interest p" implicitly in the 
operator e®^*. In order to make this dependence explicit 
one can expand e^^* = ^Lat+XQLit -^^ powers of A itera- 
tively using the relation ([5]) with A — L^ and B = XQLi: 



'o 



(17) 



e^*(Lp") = e2^*(Lp") + A / dre^^'^^^rLie^^^Lp"). +X'' fdn f 'dT2 e^°^'-^'^ QLie^°^^'-^''> QLie^°^' + 

Ja Jo Jo 
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Substitution of this into Eg. (1161) for the random force 
Kn{t) leads to the expansion 



i=0 



The lowest order term reads 

Kl^\t)=np--^F^{t), 



(18) 



(19) 



where ^o(i) = e^°*i^ is the force exerted by the bath on 
the fixed particle, which can be called the pressure force. 
As follows from and ((T7)) , the higher order terms Ki 
can be obtained recurrently as follows 

Klf>{t)^ f' dTe'^o^'-^^QLiKi^-'HT). (20) 



Note that each term in the expansion of the random force 

is zero-centered, {Kn\t)) = 0. 

In the lowest order in A the exact equation of motion 
(lisp takes the form of the generalized Langevin equation 



dt 



(21) 



+A^(n-1) / dTCo{t-T)p^-'{T), 
Jo 

where the memory kernel cq (t) is given by the correlation 
function of the pressure force. 



co{t) = n{FFo{t)). 



(22) 



Taking the average of Eq. (j21l) one obtains for the mo- 
ments 



A.(i) = (p"(t))-Pp"(i) 
the following equation 



(23) 



^A„{t) = -X'l f dTCo{t-T)A„{T) (24) 
dt m Jq 



+X'{n~l) / drco(i-r)A„_2(T). 
Jo 

In Markovian limit it takes the familiar form [l3| 

^An{t) = -X^n-foA„{t) + X^n{n - l)^7ov4„_2(i) (25) 

with the damping coefficient 

70 = — r dtc^{t) = ^ rdt{FFo{t)). (26) 
nm Jq m Jq 



Of course, Eq. ((25|) can be derived more easily from the 
phenomenological Langevin equation 



m = -X-'joP{t)+XFo{t) 



(27) 



under the assumption that Fo{t) is Gaussian noise 
An important outcome of the above microscopic deriva- 
tion is that it shows that the assumption of Gaussian 
random force is in fact unnecessary. Another advantage 
of the non-Markovian equation (|24)) is that it holds for 
any time, while its Markovian counterpart ()25p applies 
only on a time scale longer than the characteristic time 
for the decay of the correlation function co{t). 

It is easy to show that both Markovian ([25)) and non- 
Markovian ([24l) equations describe relaxation of the mo- 
ments to equilibrium values A'^ prescribed by equilibrium 
statistical mechanics, in particular 



and in general for even n 



Al^[n-l)-Al^^ 



(29) 



(equilibrium odd moments ^In+i vanish due to symme- 
try). For non-Markovian equation ()24p . thermalization 
can be proved using Laplace transformation A„(s) — 
J^dte-^'Anit), 



An{s) 



A„(0) 



+ 



s + A2(/3/m)£o(s) 
A2(n-1) 5o(s) 



(30) 



s + XHp/m)£ois) 
where Ao{s) = 1/s. For asymptotic long-time values 

= lim A„(t) = lim si„(s) (31) 

Eg. ([50)1 gives the thermal equilibrium result (1^^ . 

It also can be seen from from (|30p that relaxation to 
thermal equilibrium does not occur if 



co(s) ~ s , S>1, as s -> 0. 
In this case the damping coefficient 70 vanishes 



70 



f3 



dt co(t) 



— £o(0) =0, 



(32) 



(33) 



which corresponds to superdiffusion of the particle ■ 
Relation ((32|) as a condition of non-ergodic relaxation of 
the second moment A2{t) was discussed in [l^. As we 
see, in the lowest order in A the same condition holds for 
higher moments as well. 



IV. HIGHER-ORDER RESULTS 

As follows from (|15p . higher-order terms in the A- 
expansion of the fluctuating force 



K^{t) = K^^\t) + XKi'\t) + X' Ki^Ht) 



(34) 
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appear in the equation for the moments An in the form of 
correlations {FKn\t)). Evaluation of these correlations 
may be discouragingly complicated even for simplified 
models [l6|. However, as we show in this section, a de- 
tailed evaluation of microscopic correlations is unneces- 
sary to demonstrate convergence to thermal equilibrium 
to any perturbation order. All one actually needs is to 

find the explicit dependence of correlations {FKn \t)) on 
the particle's momentum p, which is a much easier task. 

First, recall that in the lowest order Ko{t) = 
np^-^Foit), so 

{FK'^\t))^c,{t)p-\ (35) 

with co(i) = n{FFo{t)). Substitution of this into the 
exact equation (ITSt and taking the average leads to an 
equation for the moments in the form (j24p . 

Next, it follows from the recurrence relation ([20)1 
that ~ LiKi"^ = [{m-^d/dx) p + F d/dp] K^"\ 

Therefore the correlations {FKn\t)) as functions of p 
can be obtained recurrently as follows 

{FKi^+'Ht)) r.(p + d/dp) {FKl:\t)). (36) 

From (PSl) and ([55)1 one obtains 

{FKll\t)) = c,o{t)p- + cn(i)p""'- (37) 

Explicit evaluation of the functions cio(i) and cii(t) (see 
Appendix) immediately reveals that for a homogeneous 
bath both functions vanish identically. More generally, it 
can be proved with the standard symmetry argument [17*1 
that 

{FK^^\t))^{) for odd i. (38) 

Then the first non-vanishing correction to the kernel 
{FKk°\t)) is X^{FK^n\t)). From ^ and ^ one gets 

{FKl?{t)) = C2o(i)p"+^ +C2i(t)p"-i +C22(i)p"^'- (39) 

Explicit expressions for the functions C2i(t) are not 
needed for our purposes, yet for the sake of complete- 
ness they are given in the Appendix. It is helpful, how- 
ever, to notice that the function C22{t) involves the factor 
(n — l)(n — 2), which makes it vanish for n < 3, so that 
the above expression involves only positive powers of p. 
The same is true for higher order corrections. The next 
non-zero term has the form 

{FK(^\t)) - C4o(t)p"+3 + C4i(t)p"+i (40) 

+ C42(i) + C43(t)p"-' + C44(^)p"-^ 

which can be obtained by applying twice the recurrence 
relation (|36)) to the correlation {FKn (t)) given by (p9|. 
The functions Cij(t) also depend on n. One can show 
that the above expression involves only positive powers 
of p, since C43 ~ nfc=i('^ ~ ^) ^44 ^ nfc=i("' ~ 



With the pattern given by the above relations, the gen- 
eral expression can be written in the form 

i 

{FK^^{t))^p-'Y.c.,{t)p^~^\ (41) 

where Cij (t) vanish identically for odd i due to symmetry, 
and coo{t) = co{t). We see that the higher the order of 
perturbation i, the larger the number of powers to 
which the variable of interest p" is coupled to. 

Equipped with the above relations, one can write the 
equation for the moments An{t) = (p"(t)) to any order 
in A. As an example, let us consider the equations for 
the second moment A2. In order to get the first non- 
zero correction to the lowest-order results discussed in 
the previous section, we need to expand the fiuctuating 
force K2{t) up to order A^, 

K2{t) = K^°\t) + \K^^\t) + X' Kf\t). (42) 

Then from ^ and ^ with we get 

{FK2{t))=co{t)p + \^{c2o{t)p^ + C2i{t)p]. (43) 

Substituting this into the exact equation (jlSp one obtains 
the following equation 

j^A2{t) = A^ [1 - (/3/m) A2] o Co (44) 
+A4{[3c2o-(/3/m) C21] o A2- {(3/m) C20 o A4+ C21 o 1} . 

Here and below we adopt the shorthand notation fog for 
the convolution integral dr f{T)g{t — t). One observes 
that to the given order A2 is coupled to A^, which is in 
contrast with the lowest-order approximation where the 
equation for A2 is closed. 

Applying Laplace transformation one can write the 
long-time stationary value for A2 as a fraction 

N 

lim A2{t) = lim sA2{s) = — (45) 

with the denominator 

D = X\li/m) 2o(0) + XHfi/Tn) £21(0) - iX" £20(0) 
and the numerator 

N = A2So(0) + A*C2i(0) - A^ (/3/m) £20(0) lim sAi{s). 

The stationary value for the fourth moment 
lims_).o 5^44(5) — limf_i.oo ^4(i) appears here multi- 
plied by A**, and therefore the value 3(to//3)^ found in 
the lowest-order limit, Eq. p8|) . should be assigned to 
it. Then (|45|) gives for the second moment the same 
equilibrium value as in the lowest perturbation order, 
A2{t) Al=m/p. 

Higher moments can be handled in a similar way. For 
instance, as follows from (PSj) . (1551) and ([M)) . for the 
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fourth moment A4 = (p^) the kernel {FKiit)) to order 
takes the form 

(Fi^4W) -Cop3 + A2{c2oP'Vc2lp3+C22p}. (46) 

Substitution of this into ([T5|) gives the equation 

j^M{t) = A2 [3co o A2 - {p/m) Co o Ai] (47) 

+A4{[3c2i - (/3/m) C22] o A2 + [5c2o - (P/m) C21] o A4 

-{Pirn) C20 0^6 + C22 o l|. 

Applying Laplace transform and recalling that in the low- 
est order limg^o sAq{s) = Aq = lb {m/P)^, one finds 

lim Ai{t) = lim sAi{s) = ^4 = 3 (m/pf, (48) 

which is again the equilibrium result which we already 
obtained in the lowest order. 

No new features appear as one extends the technique 
to higher perturbation orders. The next non-zero cor- 
rection corresponds to the expansion of the fluctuating 
force to order A'', Kn{t) = Ei=o^*-^"HO- The corre- 
lations {FKn\t)) for i — 0,2,4 are given by equations 
([55]) . (15^ . and (Uni), respectively. For example, for the 
second moment A2 one obtains 

{FK2{t)) = coP + A2{c2oP^ + C2ip} (49) 

-t- A'' {C40P^ -I- C4lP^ -I- C42P} • 

Then substitution into (jl5|) leads to an equation which 
differs from Eg. dl^ by the presence of terms of order A^, 

^A2{t)^\^l-{p/m)A2]oco (50) 
dt 

+A^{ [3 C20 - (/3/m) C21] 0A2- (/3/m) C20 o A4 + C21 o 1 1 
+A'^{ [3 C41 - {p/m) C42] o A2 -I- [5 C40 - iP/m) C41] o A^ 
-{P/m) C40 o Ag+ C42 o l|. 

Applying Laplace transformation and assigning equilib- 
rium values found to lower perturbation orders for long- 
time limits of A4 and Aq, one again obtains A2{t) — s- 
Al = m/p. 

V. CONCLUSION 

In the weak coupling limit, the equations for the first 
two moments Ai = (p) and A2 — (p^) of the Brownian 
particle's momentum are closed, while higher moments 
An are coupled to An~2 only, see Eq. (^5]) . To higher 
orders in the weak-coupling parameter A, a larger num- 
ber of moments are coupled. The higher perturbation 
order, the larger the number of different moments ap- 
pear in the equation for An- In a Markovian limit the 



equations coincide with those obtained by van Kampen 
within the instantaneous binary collision model, but con- 
tain parameters expressed in a totally microscopic way. 
For example, according to (HU, to order A^ the equation 
for A2 is not closed but involves coupling to A4: 

^Mt) = -liMt) - 72^4 + 73- (51) 

Dissipative coefficients 7^ are given by fluctuation- 
dissipation relations 

71 = A^(/3/m) ao - 3A^ ^20 + A^(/3/m) Q!2i, 

72 = XHP/m)a2o, (52) 

73 = ao + A^ a2i. 

Here, coefficients ao = co{t) dt and a2i = /g^ C2i{t) dt 
are system-dependent parameters, given by integrals of 
microscopic correlations. 

Since more than one microscopic parameters a are in- 
volved in (j52p . it appears natural to ask whether any 
constraints on their relations do exist which ensure relax- 
ation of the system to thermal equilibrium with the bath. 
The present paper shows that the system's thermaliza- 
tion is guaranteed by fluctuation-dissipation relations 
alone and no additional relations between microscopic 
parameters are required. Convergence to thermal equi- 
librium with the bath occurs to any order in A, in both 
Markovian and non-Markovian regimes. For instance, 
given the asymptotic result A4{t) A\ — 2>{m/jP)'^ 
found in the weak coupling limit, Eqs. ([5T|) and ([5^ 
give for A2 in the long-time limit the equilibrium value 
A2{t) ^ (73 - 72 Al)hi = m/p. 

That thermalization puts no constraints on micro- 
scopic parameters a opens an attractive avenue for phe- 
nomenological modeling. One cannot use, say, Eq. (jSip 
with arbitrary postulated values for coefficients 7^ since 
such an equation in general would disagree with equi- 
librium statistics. On the other hand, Eq.()51|) supple- 
mented with fluctuation-dissipation relations (|52p for 7^ 
with arbitrary a is thermodynamically consistent. 

Although our attention here was focused on the issue 
of thermalization and consistency with equilibrium sta- 
tistical mechanics, most interesting applications of the 
developed formalism are expected, of course, for time de- 
pendent phenomena. The coupling of a larger number of 
moments may result in much richer dynamics compared 
to that in the weak coupling limit. For a Markovian limit 
this was illustrated in [ll-Q , but results from these studies 
may be obtained (and in fact, most of them were) within 
the framework of more simple van Kampen theory. For 
future studies, it would be interesting to identify situa- 
tions where the non-Markovian form of the equations ob- 
tained in this paper would be essential and responsible for 
qualitatively new features. Application to Kramers' acti- 
vated escape problem seems particularly promising, con- 
sidering the recent demonstration that non-linear correc- 
tions to the dissipative force, which are of higher orders 
in A, may be important in the underdamped regime 
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There are several limitations of the presented study. 
One inevitable loophole is the tacit assumption that each 
term in the A-expansion is bounded for all time. To the 
best of our knowledge a general proof of this is still lack- 
ing. Also, we have assumed that the relative smallness of 
terms is determined solely by their dependence on A. For 
instance, in ([52l) 72 ^ A** is assumed to be smaller than 



71 A^ . This is not necessarily true since 71 involves fac- 
tors which are integrated correlation functions and may 
vanish identically or be very small. In such situations 
the system may exhibit non-ergodic behavior (does not 
thermilize to the bath temperature) as discussed else- 
where [l5l |. A more exotic condition of crgodicity break- 
ing is considered in [l8l |. 
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Appendix 

Exact explicit expressions for the functions Cy(t) in Eqs. (|37| and ([39|) can be written respectively as 

/ dTCij{t,T), C2j{t) ^ / dTi / dT2C2j{t,Ti,T2) 

Jo Jo Jo 



where correlation functions Cij read 

m 1^ \ ox I \ ox 

Cn = (n - 1) {{FFoit - t)Fo(<)) - (F)(FFo(t))} . 

C20 = /Fe^"(*--)^e^o(---)|-^^o(r2)\ , 
\ ox ox I 

C21 = ^ /Fe^"(*-^^)^[Fo(Ti - T2)Fo{t^)]\ + £ IfFoH - r,)e'^o{t-r,) ^p^^^^^ ^ + ^{FFo{t - n)) {FFo{t2)) 



C2 



[n - l)(n - 2)[{FFo{t - Ti)Fo{t ~ T2)Fo{t)) - {FFoit - n)){FFo{r2))} . 



Here F — F{0) — ^o(O). For a homogeneous bath Cio — Cu = by symmetry. The identity (|14l) was used for the 
derivation of the last term of C21 . Explicit evaluation of these and similar functions for a particular model is discussed 
in lia]. 
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